Off-lattice Monte Carlo Simulation of Supramolecular Polymer Architectures 



O 

<N 

o 

Q 

<N 



O 

a 

a 
o 
o 



> 

ON 
00 

m 
in 



H.E. Amuasi 1 , C. Storm 1 ' 2 
1 Department of Applied Physics and institute for Complex Molecular Systems, 
Eindhoven University of Technology, P. 0. Box 513, NL-5600 MB Eindhoven, The Netherlands 

We introduce an efficient, scalable Monte Carlo algorithm to simulate cross- linked architectures 
of freely-jointed and discrete worm-like chains. Bond movement is based on the discrete tractrix 
construction, which effects conformational changes that exactly preserve fixed-length constraints of 
all bonds. The algorithm reproduces known end-to-end distance distributions for simple, analytically 
tractable systems of cross-linked stiff and freely jointed polymers flawlessly, and is used to determine 
the effective persistence length of short bundles of semi-flexible worm-like chains, cross-linked to each 
other. It reveals a possible regulatory mechanism in bundled networks: the effective persistence of 
bundles is controlled by the linker density. 

PACS numbers: 87.16.Ka, 87.15.La, 87.16.af 



The Kratky-Porod worm- like chain (WLC) fl has 
proven to be an indispensable model for the coarse- 
grained description of stiff polymers. Biophysicists, in 
particular, have applied the model to glean the mechan- 
ics of a large variety of biological filaments, including 
actin EL Q , double-stranded DNA (3, [H| , unstructured 
RNA [6j, fibrin 0, tropocollagen U and many other 
polypeptides. However, in biologically relevant settings 
such filaments hardly ever occur or function alone as sin- 
gle chains. Instead, supramolecular associations — co- 
valent or transient — locking many chains into bundles 
or networks are the prevalent motif. While it is per- 
haps only natural to consider supramolecular structures 
of WLC's and study the effect of cross- linking in them, 
attempts to do so are severely hampered by the practical 
incompatibility of the connectivity constraints inherent 
to the architecture with adaptable and fast numerical 
algorithms. This Letter outlines an effective, scalable 
method for numerical analysis of such architectures, and 
clears the way for precise statistical-mechanical analysis 
of realistic supramolecular polymeric assemblies. 

The WLC is specified by a continuously differentiable 
space curve r(s) of length £ c parametrized by the arc 
length parameter s. It is further endowed with a Hamil- 
tonian that quantifies the cost of bending the curve: 



ds 



ds 2 



(i) 



where n is the bending modulus. Implicit in this def- 
inition is the constraint of local inextensibility, that is, 
the local tangent magnitude \dr/ds\ is unity. The per- 
sistence length £ p = n/k^T is the characteristic length 
governing the decay of tangent-tangent correlations and 
provides a quantitative measure for a polymer's flexibil- 
ity. Though the specification of the WLC model appears 
to be simple, the constraint of local inextensibility in- 
herent in the model leads to considerable mathematical 
difficulty when attempting to obtain an analytical solu- 
tion of even the simplest of such thermally fluctuating 



network structures. Nonetheless, in the so-called semi- 
flexible (£ p > £ c ) limit the radial distribution function 
may be obtained analytically @. Furthermore, isotropic 
random networks of WLC's serve as a model for the me- 
chanics of general filamentous biomaterials [3|, |7j . 

Laboratory experiments and computer simulations 
may have to pave the way to investigate the properties 
of those biological networks that remain analytically in- 
tractable. Even so, non-trivial complications arise since 
one often, as a first step, needs to discretize the WLC 
reducing it to a chain of tethers of fixed length (see, e.g., 
[Tot). The widely used Molecular Dynamics (MD) con- 
straint algorithms, such as SHAKE [Tl| and RATTLE 
[l2| have been developed to deal with these fixed-length 
constraints, but have, in practice, been limited to tree- 
structures and rigid loops — see Fig. [Ub) and c). 

Markov Chain Monte Carlo (MCMC) constraint algo- 
rithms offer a tantalizing alternative in that being purely 
stochastic, they allow for unphysical moves thus elim- 
inating the need for time-step integration and quickly 
providing good equilibrium statistics. Over the past 30 
years, a number of 'smart' MCMC moves have been ad- 
vanced for the simulation of atomistic models of melts 



of polymeric systems [13l-ll5|. After a few modifications 
(for example, after removing the fixed bond-angle con- 
straint inherent in many atomic-scale models) most of 
these techniques can be carried over into the simulation of 
the more coarse-grained discrete WLC models. However, 
these techniques are limited in the variety of cross- linked 
architectures which they are able to address 
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In this Letter we introduce TRACTRIX: a MCMC 
move that may be used to accurately simulate various 
cross-linked freely- jointed and discrete WLC architec- 
tures (Fig.HJ. The common feature of the structures (d), 
(e) and (f) in Fig. [T] that distinguishes them from the rest 
[(a), (b) and (c)] is the existence of conjoined closed loops 
that share one or more polymer links, a feature we wish 
to address, and which is a central characteristic of both 
supramolecular filaments and cross-linked networks. 

Three main technical problems are to be addressed in 




FIG. 1. (a) Linear chain; (b) loopless branched polymer; (c) 
closed-loop; (d), (e) and (f) are network structures which pos- 
sess two or more closed-loops that share the same polymeric 
link. The MCMC move TRACTRIX, introduced in this Let- 
ter, can address these structures. 



the implementation of our method: (i) preservation of 
the connectivity of the network structure, (ii) confor- 
mance to the fixed-length constraints of the inter-linking 
tethers, and (in) detailed balance. For the latter require- 
ment, we adhere to the recipe of the standard Metropo- 
lis algorithm, which is to ensure that each trial move 
is reversible, and that its probability of acceptance - 
the so-called acceptance ratio — is correctly computed 
to eventually yield Boltzmann statistics. As pointed out, 
for instance, by Maggs [ItJ, in continuum systems such 
as ours the volume element in the vicinity of the state 
is also transformed, and we must therefore consider the 
Jacobian determinant of the transformation when deter- 
mining the acceptance ratio. 

The basic unit of any network is the star which con- 
sists of n s linear chains terminating at a single central 
node (see Fig. [2]) by means of a cross-link. Any positive 
integer value for n s constitutes a star, but typically for 
biological networks n s = 4. Without loss of generality, 
let us consider the 3-arm star illustrated in Fig. EJd and 
assume that the ends A, £>, and C of its arms are tem- 
porarily fixed in space but that the central node O is 
free to move. To preserve the network connectivity at all 
times we must, whenever O is displaced by some vector 
A, move all the ends of the linear chains that terminate 
at O by the same displacement. Typically A is a ran- 
dom displacement chosen from a spherically symmetric 
distribution during the simulation. We may thus treat 
each of the linear chains independently so our focus may 
now narrow down to a single linear chain anchored at one 
end but with the other end free to move. The problem at 
hand may be set forth in two parts: firstly, given the ini- 
tial contour of the chain, how may we reversibly deform 
it so that its free end is displaced by exactly A? [18[ In 
effect, we seek a suitable invertible transform Ga that 
will act on the chain incrementing its end-to-end vector 
by A. Secondly, with what probability must we accept 
this deformation? 

We will proceed by considering a linear chain for which 
both ends are initially free to move and we will determine 
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FIG. 2. (a) Section of a network of WLC's. (b) The basic unit 
of a network: the star. When the central node O is displaced 
by a random vector A, each of the arms [e.g., link AO in (c)] 
of the star may be treated independently. 

its new position after we displace one end by So. The 
chain may be discretized into a series of TV bonds ti, t2, 
. . ., tjv of fixed length t interlinking the sequence of coor- 
dinates ro, ri, . . ., tat, so that the i-th bond ti = r^— r^_i, 
and \ri — r^_i | — t = (t = const, and i = 0, 1, . . . , N). 
The bending energy of the discretized WLC is given by 
@£({rj) = -e^iJ^U ■ t i+ i, where e = k B T£ p /t 3 . 
This expression can be shown to approach the energy for 
the continuous chain in Eq. (pQ) in the limit of N — >> oo 
and t — )• while keeping constant Nt = £ c and et 2 /N . 
However, it is sufficient to discretize the WLC so that 
there are at least 3 bonds in one persistence length 
Let gdi-i be an operator that will transform the pair 
{r^_i, r^} into {r-_ 1? r-} so that r f i _ 1 = r^_i +<5i_i and 
|r-_ x — r-| = |r^_i — r^| = t. One may readily select any 
known transformation that satisfies these two latter re- 
quirements, but for reasons that we will justify shortly, 
let us adopt Hoffman's [19| discrete tractrix construction 
[20| which is illustrated geometrically in Fig. [3£i: first, 
rigidly translate the bond ti by Si-\ so that {r^_i, r^} 
is moved to {r^_ l5 r"}. Next, form the parallelogram 
spanned by Si-\ and /t^, where / is some adjustable 
factor (/ = 2 in Fig. [3]). Finally, to obtain r[ reflect y'( in 
the parallelogram's diagonal which passes through r f i _ 1 . 
Fig. [3}d shows that following the same procedure for the 
inverse transform g^ 1 1 = g-^_ 1? but this time starting 
from the initial position {r-_ 1? r-} we recover the orig- 
inal pair {r^_i, r^} thus demonstrating the transform's 
reversibility. It can be shown that 

r- = g T (r i? ri_i, /) = p; + 2w^ ^ W ' (2) 

where pi = r-_ x - U and w; = fn + (1 - /)r»_i - r-_ x . 

We can now effect a transformation, denoted by G^ o , 
on the entire chain by displacing the polymer end ro 
by #o so that r = ro + £o> then successively apply- 
ing the discrete tractrix transformation with constant 
/ to every bond ti, t2, . . ., tjy in that order, each 
time setting Si-\ = r f i _ 1 — r^_i. Eventually, we obtain 
a new conformation of the chain: {r , r^, r^} = 
G^ o {ro, ri, . . . , tn}. Moreover we find that G^ o is re- 
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do not contribute to the value of the determinant. The 
key to computing Jq a lies in first differentiating Eq. (j3j) 

(s) 

with respect to and solving it to obtain: 
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FIG. 3. (a) Discrete tractrix transformation of the bond t$: 
1. Rigidly translate the bond by <5i_i so that {r^-i, r^} is 
moved to {r' i _ 1 , v'[ }. 2. Construct the parallelogram spanned 
by di-i and 2t». 3. Reflect r" in the parallelogram's diagonal 
(dashed line) to obtain v\. (b) Inverse of the transformation 
in (a) showing that the tractrix transformation is reversible. 
In both cases the length of the bond is preserved. 



versible since to recover the original configuration all we 
need to do is apply to the new configuration. An- 

other important feature of this transformation is that in 
general the end-to-end vector R = ro — r/v changes, al- 
beit not by So, since the other end is also displaced 
in the process. However, suppose we wanted to change 
the end-to-end vector of the polymer by exactly A, there 
may exist some So = S* for which this is possible, i.e., 
ro + A - r N = r + S* - r' N . 

Thus we may address the aforementioned problem of 
the anchored polymer chain by temporarily setting it 
free, then determining G^* for which R is incremented 
by exactly A. After applying G^*, we rigidly trans- 
late the entire chain by a displacement — r' N so 
that the tail end which was temporarily set free coin- 
cides once more with its previous position. We summa- 
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rize the full transformation {r s , r-j_ , r^_ 1 } 

{i*q S+1 \ r^ s+1 \ r^p} of the anchored linear dis- 
cretized chain as follows: 



r («+i) 



A o 

rf + (r 



5* + (r< a) 



N 



N ) 



(s) 
N 



N '), 0<i<N, 
where S* solves the second equality in Eq. (|3]), and 



A, (3) 

(4) 



/(*) 
o 



rf > = gj[ r W r f\, rf\ - r<f\, /), < i < N. (5) 



In n dimensions, Eq. (j3j) is a system of n nonlinear equa- 
tions in n unknowns which may be solved for S* . The 
acceptance ratio of G a is given by 
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the last factor being the Jacobian determinant det(Jc A ) 
of the transformation. Notice that we have used the fact 

<9r (s+1) 

= ^ofcln, where l n is the n x n identity 



that 



Or 



matrix and Sij is the kronecker-delta [see Eq. (|3J)], to 



where both derivatives in the right-hand-side may be 
found by differentiating Eq. (j2j) appropriately. The other 



derivatives 



dr 



for i = 1 , ... N follow recursively 



from the latter equation after differentiating Eq. (J5j) . Fi- 
nally, one may obtain all the matrix elements of Jq a by 
differentiating Eq. ^ and substituting these results. 

The determinant itself may be computed numerically 
by using LL^-decomposition which has a complexity of 
0(N 3 ). The algorithm is fully scalable — a suitable cut- 
off for the number of bonds TV taking part in the discrete 
tractrix move may be chosen beforehand to suit the speed 
of the computer. 

Typically, for one simulation step, a central node is 
picked at random and displaced by A. The correspond- 
ing 5* for each arm originating from the central node is 
numerically solved to a specified precision and Ga ap- 
plied. If no solution for (5* is found for an arm, or if the 
numerical solver cannot reproduce the original configu- 
ration after the inverse transformation G_a is applied, 
then the entire simulation step is rejected in accordance 
with the Metropolis algorithm rules [Hj, otherwise the 
complete acceptance ratio for the star's deformation is 
found by computing its total change in energy, the prod- 
uct of its arms' Jacobian determinants, and finally plug- 
ging them into Eq. (j6j) . To ensure ergodicity, a few steps 
with random crank-shaft rotations can be applied to each 
arm between tractrix moves — the crank-shaft rotations 
will enable each linking linear chain between nodes of the 
network to explore more of its possible conformations as 
the nodes remain fixed in space, while the discrete trac- 
trix moves will displace the nodes themselves. We will 
outline several efficient variations of TRACTRIX in a 
future publication. 

Our algorithm is the central result of this paper. To 
test it, we have applied it to various freely- jointed archi- 
tectures (e = 0) with both large and small numbers of 
bonds. In each case, TRACTRIX reproduced the equi- 
librium end-to-end distance distributions in exact agree- 
ment with their predicted analytical results. For systems 
where no such results exist, we conclude with a simple 
demonstration of TRACTRIX' proper functioning and 
use: Fig. |4] shows the results of simulations of various 
cross-linked WLC architectures, demonstrating an effect 
of some biological significance: cross-linked supramolec- 
ular polymers show a rising effective persistence length 



with cross-linking density 



Here it is defined 



eliminate the first n rows and columns of Jq a as they as the persistence length of that WLC which has the 
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FIG. 4. Simulation results (after ~ 10 7 Monte-Carlo steps): probability density distribution functions P(r) for the end-to-end 
distance r (normalized to unity) of freely fluctuating bundles of three WLC's each with £ p = 0.5£ c except for the last bundle 
which has £ p = 2£ c . Each bundle is cross- linked at equal intervals along the contour lengths of its chains. n x is the number 
of cross- links in a bundle. The red dashed curve is an analytical estimate for the last bundle: A/"47rr 2 [G f (r)] 3 where G(r) is the 
spherically symmetric radial distribution function for a stiff WLC [9| (here £ p = 2£ c ), and AT is a normalization constant. As 
expected, this prediction agrees fairly well with results for bundles with £ p > £ c . The inset shows the dependence of effective 
persistence length £ p (of those bundles with £ p = 0.5£ c ) on n x . (See text for details.) 



same £ c and expected value of the end-to-end distance 
(r) = [JdrrP(r)/(47rr 2 )]/[/drP(r)/(47rr 2 )] as the n x - 
cross-linked bundle [P(r) denotes the end-to-end distance 
distribution]. Using this design motif, Nature may cre- 
ate supramolecular filaments of tunable effective stiffness 
with only two kinds of molecules at its disposal: identi- 
cal chains and cross-linkers applied in varying concentra- 
tions. Note, too, that P(r) of an n x -cross-linked bundle 
is not the same as that of a single WLC with £ p = £p(n x ) 
— in fact, one may have to use an extensible WLC vari- 
ant to capture the complete effective mechanics. Though 
the stiffness of a bundle in reality also depends on the 
cross-linkers' stiffness and size [2l|, we did not consider 
this dependence in this initial survey. Nor did we con- 
sider the effect of excluded volume interactions between 
chains, which would result in further stiffening each bun- 
dle while incurring the additional computational cost of 
having to reject all MCMC moves that cause filaments of 
now finite cross-section to overlap or violate topological 
constraints. Extensive simulations of the collagen fibril, 
a supramolecular assembly of polypeptide triple helices 
are underway, and will be reported on in an upcoming 
publication. We emphasize that, although this work was 
inspired by bio-polymeric structures, TRACTRIX is in 
fact capable of dealing with similar configuration-space 
constraints in much more general settings and as such 
may find use well beyond biological polymers. 
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